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Abstract- In  this  study  we  compare  two  different  linear  inverse 
approaches  to  the  inverse  problem  of  electrocardiography. 
Method  A  is  the  standard  Tikhonov  zero  order  approach  used  as 
a  reference  method.  Method  B  is  based  on  the  general  deconvo¬ 
lution  theorem  proposed  by  Greensite.  Both  methods  are  applied 
to  both  the  epicardial  potential  and  the  transmembrane  poten¬ 
tial  considered  as  the  primary  source.  In  order  to  compare  the 
performance  of  both  methods  the  activation  time  pattern  is  es¬ 
timated. 
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brane  potential,  epicardial  potential,  ill-posed  problems 

I.  Introduction 

Noninvasive  functional  imaging  of  the  electrical  sources 
of  the  human  heart  is  enabled  by  electrocardiographic  (ECG) 
mapping  in  combination  with  magnetic  resonance  imaging 
(MRI).  According  to  the  bidomain  theory  [8]  the  transmem¬ 
brane  potential  cpm  is  considered  the  primary  electrical  source 
in  the  human  heart.  Due  to  admissible  assumption  of  electri¬ 
cal  isotropy  of  the  myocardium  during  the  depolarization  [2, 
8]  the  corresponding  forward  problem  degenerates  to  a  two- 
dimensional  field  and  scalar  potential  problem.  The  connec¬ 
tion  between  the  transmembrane  potential  and  the  potential 
on  all  conductivity  interfaces  as  well  as  the  body  surface  po¬ 
tential  (BSP)  is  given  by  a  Fredholm  integral  equation  of  sec¬ 
ond  kind  (Geselowitz  equation).  In  general,  the  boundary 
element  method  [1,  2]  (BEM)  is  employed  for  solving  this 
kind  of  problem  numerically  which  yields  to  an  ill-posed  and 
rank-deficient  system  of  linear  equations.  For  the  inverse 
problem  regularization  techniques  are  required  to  estimate 
parameters  describing  features  of  (pm  (e.g.  the  activation  time 
(AT)  defined  as  the  onset  of  q>m)  or  the  epicardial  potential  (pe. 
An  empirical  way  to  estimate  the  AT  from  the  epicardial  po¬ 
tential  is  presented  in  [6].  The  standard  regularization  method 
is  referred  to  as  Tikhonov  zero  order  regularization  (method 
A  in  this  study).  A  different  and  novel  method  used  in  this 
study  was  proposed  by  Greensite  [3]  (method  B). 

AT  imaging  methods  are  of  great  clinical  interest.  These 
methods  enable  the  spatially  resolved  reconstruction  of  single 
focal,  multiple  focal  and  more  complicated  activation  pat¬ 
terns.  Non-invasive  imaging  of  ectopic  and  pre-excited  ven¬ 
tricular  activation  may  be  one  of  the  possible  clinical  applica¬ 
tions. 


II.  Mathematical  Methods 

In  order  to  discretize  the  Geselowitz  equation  in  the  vol¬ 
ume  bounded  by  the  heart  surface  and  the  torso  surface  the 
BEM  was  applied  to  result  in  the  following  matrix  equation, 
when  the  transmembrane  potential  is  considered  the  source: 

D  =  L  m<t>m  (1) 


where  D  is  the  ECG  data  matrix  with  column  vectors  d,  rep¬ 
resenting  an  ECG  map  for  one  certain  time  instant  t ,  Lw  is  the 
lead-field  or  transfer  matrix  for  the  transmembrane  potential 
and  is  the  discrete  transmembrane  potential.  The  column 
vectors  of  d>w,  denoted  by  (pw,  are  the  transmembrane  poten¬ 
tial  distribution  on  the  heart  surface  at  the  time  instant  corre¬ 
sponding  to  the  column  index. 

When  the  epicardial  potential  is  the  favoured  quantity  to  re¬ 
construct  one  has  in  principal  to  solve  Laplace’s  equation  in 
the  volume  between  the  epicardial  surface  and  the  torso  sur¬ 
face.  Applying  the  BEM  yields  to  a  similar  equation  for  the 
connection  between  the  epicardial  potential  and  the  BSP: 

D  =  Le0>e  (2) 

where  Le  is  the  lead-field  matrix  for  the  discrete  epicardial 
potential  <be. 

Because  of  the  ill-posed  nature  of  the  inverse  problem,  we 
cannot  simply  invert  (1)  or  (2)  to  compute  (pw  or  (pe  Regulari¬ 
zation  schemes  which  incorporate  further  temporal  and  spa¬ 
tial  information  are  necessary.  In  the  following  section  of 
regularization  techniques  we  do  not  distinguish  formally  be¬ 
tween  the  epicardial  and  transmembrane  potential  settings.  So 
L  denotes  the  lead-field  matrix  an  cp  the  corresponding 
epicardial  or  transmembrane  potential. 


Method  A:  Standard  zero  order  Tikhonov  regularization 

For  each  fixed  value  of  t  in  T  (T  means  the  time  interval 
during  the  QRS-complex  over  which  the  ECG-data  are  col¬ 
lected),  (1)  and  (2)  reduce  to 

d,=L  <p,  (3) 

We  obtain  a  regularized  solution  to  (3)  by  minimizing  the 
following  objective  function: 

||M>«  -d,| \l  +  A0IHI2  ^min  (4) 

with  the  so-called  regularization  parameter  2.  The  first  term  in 
(4)  represents  the  least-square  solution  of  (3)  while  the  sec¬ 
ond  term  is  the  regularization  term  that  constrains  the  norm  of 
the  solution  (in  case  of  Tikhonov  zero  order).  Performing  the 
singular  value  decomposition  (S VD)  of  L, 

l  =  urv7 ,  (5) 


we  can  write  the  solution  of  (3)  as 
A(t)2 


^ =x 


r  Mt)  +<7,- 


(«,rd,k 


(6) 


were  uz  and  v*  are  the  column  vectors  of  the  matrices  U  and  V 
and  Oi  are  the  singular  values  (the  main-diagonal  elements  of 
£).  A  crucial  point  is  the  choice  of  the  optimal  regularization 
parameter  2  which  controls  the  degree  of  smoothing.  Several 
methods  for  regularization  parameter  choice  are  established 
in  literature.  In  this  study  2  is  determined  by  the  L-curve  cri¬ 
terion  [5].  This  method  involves  a  parametric  log-log  scale 
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plot  of  the  solution  norm  ||<p,||  on  the  ordinate  against  the  re¬ 
sidual  norm  ||Lcp,  -dt\\  on  the  abscissa  with  X  as  parameter. 
This  so-called  L-curve  is  plotted  for  each  point  of  time.  The 
optimal  X  is  characterized  by  the  comer  of  the  L-curve,  i.e. 
the  point  of  the  maximum  curvature.  We  emphasize  that  X  is 
computed  for  each  time  instant  separately,  hence  it  becomes 
time  dependent. 

Method  B:  improved  method 

Greensite  [3]  proposed  a  new  method  for  regularizing  the  ill- 
posed  problem  of  computing  epicardial  potentials  from  the 
body  surface  potentials.  This  method  can  also  be  utilized  to 
compute  transmembrane  potentials.  The  main  difference  to 
ordinary  Tikhonov  scheme  is  that  this  method  simultaneously 
regularizes  equations  associated  with  all  time  points.  In  the 
following  this  will  become  clearer.  Write  the  SVD  of  the  data 
matrix  as 

D  =  PQRr  ,  (7) 

where  P  is  the  matrix  of  basis  vectors  according  to  the  spatial 
domain  of  the  data,  Q  is  a  diagonal  matrix  with  the  singular 
values  as  its  diagonal  elements  and  R  is  the  matrix  of  basis 
vectors  spanning  the  time  domain  of  the  data.  The  original 
equation  transforms  now  to 

L€>  =  PQRr .  (8) 

Multiplying  (8)  by  R  and  using  the  relation  RTR=1  yields  to 

L€>R  =  PQ  .  (9) 

Let  us  denote  the  matrix  €>R  by  X  and  write  the  above  equa¬ 
tion  for  each  column  of  X  and  PQ: 

Lx,  =  <7,p; .  (10) 

Recall  that  the  pz  are  the  basis  vectors  according  to  the  spatial 
domain  of  the  data.  In  [3]  it  is  proved  that  more  accurate  re¬ 
sults  are  achieved  by  regularizing  (10),  i  .e.  to  find  the  inverse 
solution  for  the  xz  which  represent  the  sources  for  the  pz.  Us¬ 
ing  the  Tikhonov  scheme  we  obtain  the  inverse  solution 


We  reconstruct  only  equations  of  (10)  satisfying  the  discrete 
Picard  condition  [4]  while  the  other  remaining  equations  are 
eliminated  (interequation  truncation).  The  regularization  pa¬ 
rameter  is  again  determined  by  the  L-curve  criterion.  By  mul¬ 
tiplying  the  reconstructed  matrix  X  by  Rr  we  obtain  the  de¬ 
sired  potential: 

®  =  XRr  (12) 

Activation  Time 

In  order  to  compare  the  performance  of  both  methods  we 
compute  the  AT  on  the  epicardial  surface.  For  the  transmem¬ 
brane  potential  it  is  obvious  to  define  the  AT  as  time  of  the 
first  time  derivative’s  maximum  of  the  reconstructed  trans¬ 
membrane  potential  evolution.  In  case  of  the  epicardial  poten¬ 
tial  we  refer  to  the  empirical  method  presented  in  [6].  Therein 
the  AT  is  obtained  by  determining  the  time  of  occurrence  of 
the  negative  peak  of  the  first  time  derivative  of  the  potential 
at  each  source  point. 


III.  Results 

Our  simulations  were  accomplished  with  synthetic  ECG-data 
with  a  realistic  heart  and  torso  model.  The  MRI  of  the  torso 
and  heart  anatomy  was  generated  under  clinical  conditions. 
To  this  end,  the  torso  of  a  21 -year  old  male  patient  was  im¬ 
aged  with  a  1.5  Tesla  MR  scanner  (Siemens  Vision  Plus).  In 
order  to  shape  the  surface  of  the  heart,  the  cardiac  muscle  was 
imaged  in  an  ECG-gated  breath-hold  oblique  imaging  mode. 
Fig.  1.  shows  the  torso  anatomy  with  two  axial  MRI  scans, 
the  modeled  chest  surface,  the  lungs’  surfaces  and  the  model 
of  the  ventricular  surface  from  a  left-lateral  oblique  view. 


Fig.  1 .  Torso  model  generated  from  3D  MRI  data  with  two  axial  MR  scans 
from  a  left-lateral  oblique  view.  Both  lungs  and  the  ventricular  surface  are 
depicted  too.  The  spheres  represent  the  electrode  positions. 

The  epicardial  potential-to-BSP  as  well  as  the  transmembrane 
potential-to-BSP  transfer  matrix  were  computed  as  solution  to 
the  “forward  problem”.  The  triangular  mesh  for  the  epicardial 
surface  consisted  of  422  nodes. 

Synthetic  BSP-data  at  42  electrode  locations  at  the  torso’s 
surface  were  calculated  utilizing  the  forward  model.  For  this 
purpose  a  reference  AT  map  was  estimated  by  a  reference 
inverse  method  [7]  applied  to  measured  BSP  data.  This  data 
was  acquired  under  clinical  conditions  and  reflect  a  cardiac 
rhythm  resulting  from  catheter  pacing  in  the  apex  of  the  right 
ventricle.  The  reference  AT  map  is  shown  in  Fig.  2.  from  an 
anterior  epicardial  view.  Gaussian  noise  was  added  to  the 
synthetic  BSP  data  for  signal-to-noise  ratio  of  30  dB.  The 
time  points  are  taken  at  1ms  time  increments  to  encompass  a 
135ms  interval  containing  the  activation  (QRS)  interval. 
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Fig.2.  Reference  AT  map  on  the  anterior  epicardium  in  ms  estimated  from 
measured  data  by  the  method  presented  in  [7].  The  pattern  shows  more  dis¬ 
tributed  events  including  one  early  breakthrough  near  the  apex  and  two  early 
activation  events  on  left  and  right  lateral  site  of  the  anterior  epicardium. 


For  inverse  method  A  the  regularization  parameter  2  was  de¬ 
termined  at  135  time  points  by  the  L-curve  criterion  (maxi¬ 
mum  curvature).  We  experienced  that  slightly  smoothing  of 
the  temporal  2-trend  by  cubic  splines  leads  to  better  results 
for  the  inverse  solution,  in  fact,  reflected  by  an  improvement 
of  the  correlation  coefficient  for  the  AT  map  from  0.77  (non 
smoothed)  up  to  0.91  (smoothed).  The  reconstructed  AT 
maps  are  depictured  in  Fig.  3.  The  upper  panel  shows  the  AT 
result  for  the  anterior  epicardium  calculated  in  the  setting  for 
the  epicardial  potential.  In  order  to  estimate  the  AT  from  the 
reconstructed  epicardial  potential  the  first  time  derivative  had 
to  be  determined.  The  numerically  computing  of  the  time 
derivative  was  enabled  by  smoothing  with  cubic  spline  inter¬ 
polation.  A  correlation  coefficient  of  0.91  and  a  relative  error 
of  0.16  could  be  attained.  (Relative  error  is  the  square-root  of 
the  sum  of  the  squares  of  the  difference  between  reference 
and  computed  AT  divided  by  the  square-root  of  the  sum  of 
the  squares  of  the  reference  AT,  where  the  sums  are  over  all 
source  points  of  the  epicardial  surface).  In  comparison  the 
lower  panel  of  Fig.  3  reflects  the  results  of  the  transmembrane 
potential  setting.  A  similar  procedure  for  AT  computing 
yielded  to  values  of  0.88  and  0.17  for  the  correlation  coeffi¬ 
cient  and  the  relative  error,  respectively.  In  view  of  both  re¬ 
sults  we  can  not  state  essential  differences  using  the  epicar¬ 
dial  potential  formalism  and  the  transmembrane  potential 
formalism.  The  main  characteristics  of  the  reference  AT  pat¬ 
tern  could  be  retrieved  in  both  approaches. 

For  inverse  method  B  equation  sequence  (10)  constitutes 
42  equations  (i=l,  ...  ,42)  according  to  42  basis  vectors  p,  of 
the  spatial  domain  on  the  electrodes  site.  As  mentioned  before 
it  is  meaningful  to  solve  the  mth  equation  only  if  the  discrete 
Picard  condition  is  satisfied.  Thus,  for  the  first  8  equations  the 


Fig.3.  Reconstructed  AT  map  on  the  anterior  epicardium  applying  inverse 
method  A.  The  upper  panel  shows  the  solution  based  on  the  epicardial  setting 
in  comparison  to  the  lower  panel,  where  the  transmembrane  potential  repre¬ 
sents  the  source. 

solution  was  obtained  by  (11)  calculating  the  suited  regulari¬ 
zation  parameter  by  the  L-curve  criterion.  The  results  for  the 
AT  are  presented  in  Fig.  4.  In  the  upper  panel  the  estimated 
AT  in  the  epicardial  potential  setting  is  illustrated  The  corre¬ 
sponding  correlation  coefficient  is  disposed  to  0.90  and  for 
the  relative  error  a  value  of  0.17  was  obtained.  In  addition, 
the  result  of  the  transmembrane  potential  formulation  is 
shown  in  the  lower  panel.  Similar  values  of  0.91  and  0.17  for 
the  correlation  coefficient  and  the  relative  error,  respectively, 
can  be  found.  Again,  no  substantial  varieties  in  both,  the 
epicardial  and  transmembrane  potential  setting  can  be  seen 
from  Fig.  4.  Both  display  the  main  events  of  the  reference 
pattern  in  an  acceptable  way.  However,  in  contrast 
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brane  potential  formalism  which  allows  to  reconstruct  source 
patterns  on  the  epicardium  and  on  the  endocardium  as  well. 
On  the  epicardium  both  settings  lead  to  similar  results  as 
shown  in  this  study.  The  transmembrane  potential  formalism, 
however,  enables  more  powerful  regularization  methods,  be¬ 
cause  more  physiological  a-priori  information  can  be  imposed 
on  the  sought  solution.  Therefore  further  research  will  be 
focused  on  the  spatio-temporal  regularization  of  the  trans¬ 
membrane  potential  inverse  approach. 
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Fig.4.  Reconstructed  AT  map  on  the  anterior  epicardium  applying  inverse 
method  B.  The  upper  panel  shows  the  solution  based  on  the  epicardial  setting 
in  comparison  to  the  lower  panel,  where  the  transmembrane  potential  repre¬ 
sents  the  source. 

to  method  A  method  B  yields  better  results,  which  can  be 
seen,  in  particular,  by  smoother  AT  patterns. 

IV.  Discussion  and  Conclusions 


The  reconstruction  of  epicardial  potentials  from  the  BSP  con¬ 
stitutes  the  standard  approach  to  the  electrocardiographic  in¬ 
verse  problem.  To  provide  the  cardiologist  with  information 
of  clinical  interest  (e.g.,  the  AT),  however,  one  has  to  solve 
another  inverse  problem  or  to  utilize  some  empirical  methods 
as  done  in  this  study.  Furthermore,  no  information  regarding 
the  endocardium  can  be  extracted  with  the  epicardial  setting. 
To  overcome  this  drawback  one  can  employ  the  transmem- 


